oneimpact pacakge# libraries
library(rgrass7)
library(raster)
library(terra)
library(NinaR)
library(dplyr)
library(viridis)
library(oneimpact)
Connect to GRASS.
# start GRASS and connect to my mapset
grassdir <- system("grass78 --config path", intern = T)
gisDB <- "/data/grass"
loc <- "ETRS_33N/"
ms <- "u_bb_cuminf"
# initGRASS(gisBase = grassdir,
# home = tempdir(),
# override = T,
# gisDbase = gisDB,
# location = loc,
# mapset = ms)
# more directly within NINA
NinaR::grassConnect(mapset = ms)
## gisdbase /data/grass
## location ETRS_33N
## mapset u_bb_cuminf
## rows 361
## columns 478
## north 6658900
## south 6622800
## west 146900
## east 194700
## nsres 100
## ewres 100
# copy test region vector and map of private cabins
# check region
gmeta() # g.copy is not affected by the region
## gisdbase /data/grass
## location ETRS_33N
## mapset u_bb_cuminf
## rows 361
## columns 478
## north 6658900
## south 6622800
## west 146900
## east 194700
## nsres 100
## ewres 100
# copy region vector
rgrass7::execGRASS("g.copy", parameters = list(vector = "region_test_influence@u_bernardo.brandao,region_test_influence"),
flags = c("overwrite"))
## Warning in rgrass7::execGRASS("g.copy", parameters = list(vector = "region_test_influence@u_bernardo.brandao,region_test_influence"), : The command:
## g.copy --overwrite vector=region_test_influence@u_bernardo.brandao,region_test_influence
## produced at least one warning during execution:
## Copy vector <region_test_influence@u_bernardo.brandao> to current mapset as
## <region_test_influence>
## WARNING: Vector map <region_test_influence> already exists and will be
## overwritten
## Copy vector <region_test_influence@u_bernardo.brandao> to current mapset as
## <region_test_influence>
## WARNING: Vector map <region_test_influence> already exists and will be
## overwritten
# copy private cabins
rgrass7::execGRASS("g.copy", parameters = list(raster = "private_cabins_rast@u_bernardo.brandao,private_cabins_rast"),
flags = c("overwrite"))
## Copy raster <private_cabins_rast@u_bernardo.brandao> to current mapset as
## <private_cabins_rast>
# define region
# I defined a region_test_influence through the GRASS GUI and using v.in.region output=region_test_influence
rgrass7::execGRASS("g.region", parameters = list(vector = "region_test_influence", align = "private_cabins_rast"),
flags = "p")
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
# create a new small map for the current region, just to ease visualization in R
rgrass7::execGRASS("r.mapcalc", expression = "private_cabins_sub = if(!isnull(private_cabins_rast), 1, null())",
flags = c("overwrite", "quiet"))
# show input map and region here
rgrass7::use_sp()
cabins <- readRAST(c("private_cabins_rast_sub")) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_rast_sub has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
terra::plot(cabins, col = "black")
# Euclidean distance
# R
cabins_dist_r <- calc_influence_nearest(cabins, where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_dist_names <- calc_influence_nearest(name_var, where = "GRASS", quiet = F, overwrite = T)
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## Writing output raster maps...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
# visualize
cabins_dist_grass <- readRAST(cabins_dist_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_dist_r, cabins_dist_grass),
main = paste("Euclidean distance from private cabins - ",
c("R", "GRASS")))
# Log distance
# R
cabins_logdist_r <- calc_influence_nearest(cabins, transform = "log", log_base = 10, where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_logdist_names <- calc_influence_nearest(name_var, transform = "log", log_base = 10,
where = "GRASS", quiet = F, overwrite = T)
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## Writing output raster maps...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## [1] "log(A + 1.000000, 10.000000)"
## Calculating log-distance...
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
# visualize
cabins_logdist_grass <- readRAST(cabins_logdist_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_log has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_logdist_r, cabins_logdist_grass),
main = paste("Log-distance from private cabins - ",
c("R", "GRASS")))
# R
cabins_sqrtdist_r <- calc_influence_nearest(cabins, transform = "sqrt", where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_sqrtdist_names <- calc_influence_nearest(name_var, transform = "sqrt",
where = "GRASS", quiet = F, overwrite = T)
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## Writing output raster maps...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## [1] "sqrt(A + 1.000000)"
## Calculating sqrt-distance...
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
# visualize
cabins_sqrtdist_grass <- readRAST(cabins_sqrtdist_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_sqrt has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_sqrtdist_r, cabins_sqrtdist_grass),
main = paste("Sqrt-distance from private cabins - ",
c("R", "GRASS")))
Example with ZoI = 1000m
# R
cabins_exp_decay_r <- calc_influence_nearest(cabins, transform = "exp_decay",
zoi = 1000, where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_exp_decay_names <- calc_influence_nearest(name_var, transform = "exp_decay",
zoi = 1000, where = "GRASS",
quiet = F, overwrite = T)
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## Writing output raster maps...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## [1] "1.000000 * exp(-0.002773 * A)"
## Calculating exponential decay influence...
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
# visualize
cabins_exp_decay_grass <- readRAST(cabins_exp_decay_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_exp_decay1000 has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_exp_decay_r, cabins_exp_decay_grass),
main = paste("Exp-decay influence from private cabins - ",
c("R", "GRASS")))
Example with ZoI = 1000m
# R
cabins_bartlett_r <- calc_influence_nearest(cabins, transform = "bartlett",
zoi = 1000, where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_bartlett_names <- calc_influence_nearest(name_var, transform = "bartlett",
zoi = 1000, where = "GRASS",
quiet = F, overwrite = T)
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
## Calculating Euclidean distance...
## Reading raster map <private_cabins_sub>...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## Writing output raster maps...
## 0% 3% 6% 9% 12% 15% 18% 21% 24% 27% 30% 33% 36% 39% 42% 45% 48% 51% 54% 57% 60% 63% 66% 69% 72% 75% 78% 81% 84% 87% 90% 93% 96% 99% 100%
## [1] "if(A <= 1000.000000, 1 - (1/1000.000000) * A, 0)"
## Calculating Bartlett influence...
## projection: 1 (UTM)
## zone: 33
## datum: etrs89
## ellipsoid: grs80
## north: 6658900
## south: 6622800
## west: 146900
## east: 194700
## nsres: 100
## ewres: 100
## rows: 361
## cols: 478
## cells: 172558
# visualize
cabins_bartlett_grass <- readRAST(cabins_bartlett_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_bartlett1000 has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_bartlett_r, cabins_bartlett_grass),
main = paste("Linear-decay influence from private cabins - ",
c("R", "GRASS")))
Example with ZoI = 1000m
# R
cabins_hnorm_decay_r <- calc_influence_nearest(cabins, transform = "half_norm",
zoi = 1000, where = "R")
# GRASS
name_var <- "private_cabins_sub"
cabins_hnorm_decay_names <- calc_influence_nearest(name_var, transform = "half_norm",
zoi = 1000, where = "GRASS",
quiet = T, overwrite = T)
## Calculating Euclidean distance...
## [1] "1.000000 * exp(-0.000003 * pow(A, 2))"
## Calculating half-normal decay influence...
# visualize
cabins_exp_decay_grass <- readRAST(cabins_hnorm_decay_names) %>%
raster::raster() %>%
terra::rast()
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
## Warning in rgdal::GDALinfo(fname, silent = ignore.stderr): statistics not
## supported by this driver
## /data/grass/ETRS_33N/u_bb_cuminf/cellhd/private_cabins_sub_inf_nearest_half_norm1000 has GDAL driver GRASS
## and has 361 rows and 478 columns
## Warning in getProjectionRef(x, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = OVERRIDE_PROJ_DATUM_WITH_TOWGS84, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
## but +towgs84= values preserved
terra::plot(c(cabins_hnorm_decay_r, cabins_exp_decay_grass),
main = paste("Halfnormal influence from private cabins - ",
c("R", "GRASS")))
# Euclidean distance
# name_var <- "private_cabins_sub"
# execGRASS("r.resamp.filter", input = name_var, output = "test_neighborhood", filter = "gauss,box", radius = c(500,1500))
#
# execGRASS("r.mapcalc", expression = "private_cabins_binary = if(isnull(private_cabins_sub), 0, 1)", flags = "overwrite")
#
# execGRASS("r.neighbors", input = "private_cabins_binary", output = "test_neighborhood", size = 5, flags = "overwrite")
#
# Zone of Influence
# cumulative zone of influence
# cumulative affected area
# landcape effect, cumulative weight - "Miguet et al"
# Perceived disturbance
#
# multiply two layers referreing to two different scales of a variable density
#
# # visualize
# cabins_dens <- readRAST("test_neighborhood") %>%
# raster::raster() %>%
# terra::rast()
#
# terra::plot(cabins_dens, main = "Density")
#
# threshold = 0.05
# execGRASS("r.mapcalc", expression = "areas = if(test_neighborhood > 0.05, 1, 0)", flags = "overwrite")
#
# cabins_areas_affected <- readRAST("areas") %>%
# raster::raster() %>%
# terra::rast()
#
# terra::plot(cabins_areas_affected, main = "areas affected")
#
# #
# execGRASS("g.extension", extension = "r.area")
#
# execGRASS("r.clump", input = "areas", output = "patches", flags = "overwrite")
#
# cabins_areas_patches <- readRAST("patches") %>%
# raster::raster() %>%
# terra::rast()
#
# terra::plot(cabins_areas_affected, main = "patch areas affected")